EORTC HNSCC ICI Honore et al_02032025

library(swimplot) library(grid) library(gtable) library(readr) library(mosaic) library(dplyr) library(survival) library(survminer) library(ggplot2) library(scales) library(coxphf) library(ggthemes) library(tidyverse) library(gtsummary) library(flextable) library(parameters) library(car) library(ComplexHeatmap) library(tidyverse) library(readxl) library(janitor) library(DT) library(pROC) library(rms)

#ctDNA Detection Rates by Window and Stages

#ctDNA at Baseline
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data$ctDNA.Base <- factor(circ_data$ctDNA.Base, levels=c("NEGATIVE","POSITIVE"))
circ_data <- subset(circ_data, ctDNA.Base %in% c("NEGATIVE", "POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("O","I","II","III","IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.Base == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.Base, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
  Stage = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.Base == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Stage = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA post-treatment
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels=c("NEGATIVE","POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("O","I","II","III","IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.postTx == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.postTx, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
  Stage = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.postTx == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Stage = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#Overview plot

setwd("~/Downloads") 
clinstage <- read.csv("EORTC ICI_OP.csv")
clinstage_df <- as.data.frame(clinstage)

# Creating the basic swimmer plot
oplot <- swimmer_plot(df=clinstage_df,
                      id='PatientName',
                      end='fu.diff.months',
                      fill='gray',
                      width=.01)

# Adding themes and scales
oplot <- oplot + theme(panel.border = element_blank())
oplot <- oplot + scale_y_continuous(breaks = seq(0, 48, by = 3))
oplot <- oplot + labs(x ="Patients", y="Months from Immunotherapy Start")

# Adding swimmer points
oplot_ev1 <- oplot + swimmer_points(df_points=clinstage_df,
                                    id='PatientName',
                                    time='date.diff.months',
                                    name_shape ='Event_type',
                                    name_col = 'Event',
                                    size=3.5,fill='black')
# Optionally uncomment and use col='darkgreen' if needed

# Adding shape manual scale
oplot_ev1.1 <- oplot_ev1 + ggplot2::scale_shape_manual(name="Event_type",
                                                       values=c(1,16,6,4),
                                                       breaks=c('ctDNA_neg','ctDNA_pos','Imaging','Death'))

# Display the plot
oplot_ev1.1

                                    oplot_ev2 <- oplot_ev1.1 + swimmer_lines(df_lines=clinstage_df,
                                                                             id='PatientName',
                                                                             start='Tx_start.months',
                                                                             end='Tx_end.months',
                                                                             name_col='Tx_type',
                                                                             size=3.5,
                                                                             name_alpha = 1.0)
                                    oplot_ev2 <- oplot_ev2 + guides(linetype = guide_legend(override.aes = list(size = 5, color = "black")))
                                    oplot_ev2

oplot_ev2.2 <- oplot_ev2 + ggplot2::scale_color_manual(name="Event",values=c( "orange", "black", "black", "lightblue", "green", "red"))
oplot_ev2.2

#Overview plot - stratified by BOR

setwd("~/Downloads") 
clinstage <- read.csv("EORTC ICI_OP.csv")
clinstage_df <- as.data.frame(clinstage)

# Creating the basic swimmer plot
oplot_stratify <-swimmer_plot(df=clinstage_df,
                              id='PatientName',
                              end='fu.diff.months',
                              col="gray",
                              alpha=0.75,
                              width=.01,
                              base_size = 14,
                              stratify= c('RECIST'))
oplot_stratify <- oplot_stratify + theme(panel.border = element_blank())
oplot_stratify <- oplot_stratify + scale_y_continuous(breaks = seq(0, 42, by = 3))
oplot_stratify <- oplot_stratify + labs(x ="Patients" , y="Months from Immunotherapy Start")

# Adding swimmer points
oplot_ev1 <- oplot_stratify + swimmer_points(df_points=clinstage_df,
                                    id='PatientName',
                                    time='date.diff.months',
                                    name_shape ='Event_type',
                                    name_col = 'Event',
                                    size=3.5,fill='black')
# Optionally uncomment and use col='darkgreen' if needed

# Adding shape manual scale
oplot_ev1.1 <- oplot_ev1 + ggplot2::scale_shape_manual(name="Event_type",
                                                       values=c(1,16,6,4),
                                                       breaks=c('ctDNA_neg','ctDNA_pos','Imaging','Death'))

# Display the plot
oplot_ev1.1

                                    oplot_ev2 <- oplot_ev1.1 + swimmer_lines(df_lines=clinstage_df,
                                                                             id='PatientName',
                                                                             start='Tx_start.months',
                                                                             end='Tx_end.months',
                                                                             name_col='Tx_type',
                                                                             size=3.5,
                                                                             name_alpha = 1.0)
                                    oplot_ev2 <- oplot_ev2 + guides(linetype = guide_legend(override.aes = list(size = 5, color = "black")))
                                    oplot_ev2

oplot_ev2.2 <- oplot_ev2 + ggplot2::scale_color_manual(name="Event",values=c( "orange", "black", "black", "lightblue", "green", "red"))
oplot_ev2.2

#Association of Baseline ctDNA MTM levels with clinicopathological factors

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_datadf <- as.data.frame(circ_data)

tally(~cStage, data=circ_data, margins = TRUE)

cStage
  0-II III-IV  Total 
     7     22     29 

circ_data$cStage <- factor(circ_data$cStage, levels = c("0-II","III-IV"), labels = c("0-II (n=7)","III-IV (n=22)"))
boxplot(ctDNA.Base.MTM~cStage, data=circ_data, main="ctDNA pre-treatment MTM - Stage", xlab="Stage", ylab="MTM/mL", col="white",border="black")

m1<-wilcox.test(ctDNA.Base.MTM ~ cStage, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m1)

    Wilcoxon rank sum test with continuity correction

data:  ctDNA.Base.MTM by cStage
W = 39.5, p-value = 0.05896
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -5.19000e+01  4.27605e-05
sample estimates:
difference in location 
             -6.799975 

tally(~Inclusion.status, data=circ_data, margins = TRUE)

Inclusion.status
Loco-regional    Metastatic         Total 
           18            11            29 

circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional","Metastatic"), labels = c("Loco-regional (n=18)","Metastatic (n=11)"))
boxplot(ctDNA.Base.MTM~Inclusion.status, data=circ_data, main="ctDNA pre-treatment MTM - Disease Status", xlab="Disease Status", ylab="MTM/mL", col="white",border="black")

m2<-wilcox.test(ctDNA.Base.MTM ~ Inclusion.status, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m2)

    Wilcoxon rank sum test with continuity correction

data:  ctDNA.Base.MTM by Inclusion.status
W = 93.5, p-value = 0.8219
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -8.399949  6.899953
sample estimates:
difference in location 
            -0.3727674 

tally(~cT.Status, data=circ_data, margins = TRUE)

cT.Status
cT0-T2 cT3-T4    cTx  Total 
    11     17      1     29 

circ_data$cT.Status <- factor(circ_data$cT.Status, levels = c("cT0-T2","cT3-T4"), labels = c("cT0-T2 (n=11)","cT3-T4 (n=17)"))
boxplot(ctDNA.Base.MTM~cT.Status, data=circ_data, main="ctDNA pre-treatment MTM - cT status", xlab="cT status", ylab="MTM/mL", col="white",border="black")

m3<-wilcox.test(ctDNA.Base.MTM ~ cT.Status, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m3)

    Wilcoxon rank sum test with continuity correction

data:  ctDNA.Base.MTM by cT.Status
W = 45.5, p-value = 0.02523
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -57.0000356  -0.7999756
sample estimates:
difference in location 
             -7.400021 

tally(~cN.Status, data=circ_data, margins = TRUE)

cN.Status
   cN0 cN1-N3    cNx  Total 
    12     15      2     29 

circ_data$cN.Status <- factor(circ_data$cN.Status, levels = c("cN0","cN1-N3"), labels = c("cN0 (n=12)","cN1-N3 (n=15)"))
boxplot(ctDNA.Base.MTM~cN.Status, data=circ_data, main="ctDNA pre-treatment MTM - cN status", xlab="cN status", ylab="MTM/mL", col="white",border="black")

m4<-wilcox.test(ctDNA.Base.MTM ~ cN.Status, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m4)

    Wilcoxon rank sum test with continuity correction

data:  ctDNA.Base.MTM by cN.Status
W = 32.5, p-value = 0.005379
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -60.200006  -3.299997
sample estimates:
difference in location 
             -10.15193 

tally(~p16.status, data=circ_data, margins = TRUE)

p16.status
Negative Positive  Unknown    Total 
      24        4        1       29 

circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative","Positive"), labels = c("p16 neg (n=24)","p16 pos (n=4)"))
boxplot(ctDNA.Base.MTM~p16.status, data=circ_data, main="ctDNA pre-treatment MTM - p16 status", xlab="p16 status", ylab="MTM/mL", col="white",border="black")

m5<-wilcox.test(ctDNA.Base.MTM ~ p16.status, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m5)

    Wilcoxon rank sum test with continuity correction

data:  ctDNA.Base.MTM by p16.status
W = 55.5, p-value = 0.6453
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -47.89994  46.29998
sample estimates:
difference in location 
             0.8901786 

#Median MTM/mL levels for ctDNA positive pts pre-treatment

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]

median_ctDNA <- median(circ_data$ctDNA.Base.MTM, na.rm = TRUE)
range_ctDNA <- range(circ_data$ctDNA.Base.MTM, na.rm = TRUE)
cat("Median MTM/mL post-treatment:", median_ctDNA, "\n")

Median MTM/mL post-treatment: 8.4 

cat("Range MTM/mL post-treatment:", range_ctDNA[1], "-", range_ctDNA[2], "\n")

Range MTM/mL post-treatment: 0.2 - 986 

#Median MTM/mL levels for ctDNA positive pts post-treatment

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_data <- circ_data[circ_data$ctDNA.postTx=="POSITIVE",]

median_ctDNA <- median(circ_data$ctDNA.postTx.MTM, na.rm = TRUE)
range_ctDNA <- range(circ_data$ctDNA.postTx.MTM, na.rm = TRUE)
cat("Median MTM/mL post-treatment:", median_ctDNA, "\n")

Median MTM/mL post-treatment: 7.9 

cat("Range MTM/mL post-treatment:", range_ctDNA[1], "-", range_ctDNA[2], "\n")

Range MTM/mL post-treatment: 0.1 - 737.8 

#Median time from end treatment to progression for ctDNA negative pts post-treatment

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_data <- circ_data[circ_data$ctDNA.postTx=="NEGATIVE",]
circ_data <- circ_data[circ_data$PFS.Event=="TRUE",]

median_PFS <- median(circ_data$PFS.months, na.rm = TRUE)
range_PFS <- range(circ_data$PFS.months, na.rm = TRUE)
cat("Median PFS:", median_PFS, "\n")

Median PFS: 4.977494 

cat("Range PFS:", range_PFS[1], "-", range_PFS[2], "\n")

Range PFS: 3.055492 - 19.67999 

#PFS by ctDNA status post/during-ICI

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)~ctDNA.postTx, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) ~ 
    ctDNA.postTx, data = circ_data)

                       n events median 0.95LCL 0.95UCL
ctDNA.postTx=NEGATIVE  9      4     NA    6.70      NA
ctDNA.postTx=POSITIVE 20     19   3.81    2.37    6.64

event_summary <- circ_data %>%
  group_by(ctDNA.postTx) %>%
  summarise(
    Total = n(),
    Events = sum(PFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.postTx, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="PFS - ctDNA Post/Under treatment", ylab= "Progression-Free Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(0, 12, 24, 36))

Call: survfit(formula = surv_object ~ ctDNA.postTx, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.postTx=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0      9       0    1.000   0.000        1.000        1.000
   12      5       3    0.667   0.157        0.282        0.878
   24      3       1    0.533   0.173        0.177        0.796
   36      2       0    0.533   0.173        0.177        0.796

                ctDNA.postTx=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     20       0     1.00  0.0000      1.00000        1.000
   12      2      18     0.10  0.0671      0.01698        0.272
   24      1       1     0.05  0.0487      0.00345        0.205
   36      1       0     0.05  0.0487      0.00345        0.205

circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.postTx, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.postTx, data = circ_data)

  n= 29, number of events= 23 

                       coef exp(coef) se(coef)     z Pr(>|z|)   
ctDNA.postTxPOSITIVE 1.5616    4.7666   0.5703 2.738  0.00618 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                     exp(coef) exp(-coef) lower .95 upper .95
ctDNA.postTxPOSITIVE     4.767     0.2098     1.559     14.58

Concordance= 0.649  (se = 0.049 )
Likelihood ratio test= 9.7  on 1 df,   p=0.002
Wald test            = 7.5  on 1 df,   p=0.006
Score (logrank) test = 8.78  on 1 df,   p=0.003

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 4.77 (1.56-14.58); p = 0.006"

circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$PFS.Event <- factor(circ_data$PFS.Event, levels = c("FALSE", "TRUE"), labels = c("No Progression", "Progression"))
contingency_table <- table(circ_data$ctDNA.postTx, circ_data$PFS.Event)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 6.8323, df = 1, p-value = 0.008952

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.005482
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
    1.68957 1169.09723
sample estimates:
odds ratio 
  20.33413 

print(contingency_table)

          
           No Progression Progression
  Negative              5           4
  Positive              1          19

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA status post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "Progression",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("No Progression" = "blue", "Progression" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#OS by ctDNA status post/during-ICI

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$OS.months, event = circ_data$OS.Event)~ctDNA.postTx, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.months, event = circ_data$OS.Event) ~ 
    ctDNA.postTx, data = circ_data)

                       n events median 0.95LCL 0.95UCL
ctDNA.postTx=NEGATIVE  9      4     NA   16.85      NA
ctDNA.postTx=POSITIVE 20     17   10.4    7.75      18

event_summary <- circ_data %>%
  group_by(ctDNA.postTx) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.postTx, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA Post/Under treatment", ylab= "Overall Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(12, 24, 36))

Call: survfit(formula = surv_object ~ ctDNA.postTx, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.postTx=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      7       1    0.889   0.105        0.433        0.984
   24      3       3    0.508   0.177        0.157        0.781
   36      2       0    0.508   0.177        0.157        0.781

                ctDNA.postTx=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      9      11     0.45  0.1112       0.2311        0.647
   24      3       5     0.18  0.0900       0.0480        0.380
   36      2       1     0.12  0.0775       0.0213        0.311

circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.postTx, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.postTx, data = circ_data)

  n= 29, number of events= 21 

                       coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.postTxPOSITIVE 1.2213    3.3918   0.5621 2.173   0.0298 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                     exp(coef) exp(-coef) lower .95 upper .95
ctDNA.postTxPOSITIVE     3.392     0.2948     1.127     10.21

Concordance= 0.638  (se = 0.048 )
Likelihood ratio test= 5.81  on 1 df,   p=0.02
Wald test            = 4.72  on 1 df,   p=0.03
Score (logrank) test = 5.29  on 1 df,   p=0.02

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 3.39 (1.13-10.21); p = 0.03"

circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$OS.Event <- factor(circ_data$OS.Event, levels = c("FALSE", "TRUE"), labels = c("Alive", "Deceased"))
contingency_table <- table(circ_data$ctDNA.postTx, circ_data$OS.Event)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 3.2819, df = 1, p-value = 0.07005

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.0667
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  0.8624077 62.5051555
sample estimates:
odds ratio 
  6.505452 

print(contingency_table)

          
           Alive Deceased
  Negative     5        4
  Positive     3       17

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA status post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "Living Status",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Alive" = "blue", "Deceased" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#Association of ctDNA status post/during-ICI with BOR

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$RESIST <- factor(circ_data$RESIST, levels = c("CR", "PR", "SD", "PD"))
contingency_table <- table(circ_data$ctDNA.postTx, circ_data$RESIST)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test

data:  contingency_table
X-squared = 11.869, df = 3, p-value = 0.007847

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.006541
alternative hypothesis: two.sided

print(contingency_table)

          
           CR PR SD PD
  Negative  4  3  0  2
  Positive  1  2  7 10

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA status post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "BOR",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("CR" = "blue", "PR" = "lightblue", "SD" = "orange", "PD" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#PFS by ctDNA kinetics post/during-ICI

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ΔctDNA!="",]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)~ΔctDNA, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) ~ 
    ΔctDNA, data = circ_data)

                 n events median 0.95LCL 0.95UCL
ΔctDNA=NEGATIVE 13      8  18.04     6.7      NA
ΔctDNA=POSITIVE 12     12   2.41     2.3      NA

event_summary <- circ_data %>%
  group_by(ΔctDNA) %>%
  summarise(
    Total = n(),
    Events = sum(PFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
KM_curve <- survfit(surv_object ~ ΔctDNA, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="PFS - ctDNA Kinetics Post/Under treatment", ylab= "Progression-Free Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("Decreased ΔctDNA", "Increased ΔctDNA"), legend.title="")

summary(KM_curve, times= c(0, 12, 24, 36))

Call: survfit(formula = surv_object ~ ΔctDNA, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ΔctDNA=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     13       0    1.000   0.000        1.000        1.000
   12      6       6    0.538   0.138        0.248        0.760
   24      3       2    0.359   0.139        0.117        0.613
   36      2       0    0.359   0.139        0.117        0.613

                ΔctDNA=POSITIVE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           12            0            1            0            1            1 

circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ΔctDNA, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ΔctDNA, data = circ_data)

  n= 25, number of events= 20 

                 coef exp(coef) se(coef)    z Pr(>|z|)    
ΔctDNAPOSITIVE  2.716    15.112    0.794 3.42 0.000626 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

               exp(coef) exp(-coef) lower .95 upper .95
ΔctDNAPOSITIVE     15.11    0.06617     3.188     71.64

Concordance= 0.716  (se = 0.044 )
Likelihood ratio test= 17.9  on 1 df,   p=2e-05
Wald test            = 11.7  on 1 df,   p=6e-04
Score (logrank) test = 18.33  on 1 df,   p=2e-05

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 15.11 (3.19-71.64); p = 0.001"

circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels = c("NEGATIVE", "POSITIVE"), labels = c("Decreased", "Increased"))
circ_data$PFS.Event <- factor(circ_data$PFS.Event, levels = c("FALSE", "TRUE"), labels = c("No Progression", "Progression"))
contingency_table <- table(circ_data$ΔctDNA, circ_data$PFS.Event)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 3.6158, df = 1, p-value = 0.05723

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.03913
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.013281      Inf
sample estimates:
odds ratio 
       Inf 

print(contingency_table)

           
            No Progression Progression
  Decreased              5           8
  Increased              0          12

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA kinetics post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "Progression",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("No Progression" = "blue", "Progression" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#OS by ctDNA kinetics post/during-ICI

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ΔctDNA!="",]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$OS.months, event = circ_data$OS.Event)~ΔctDNA, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.months, event = circ_data$OS.Event) ~ 
    ΔctDNA, data = circ_data)

                 n events median 0.95LCL 0.95UCL
ΔctDNA=NEGATIVE 13      8   18.0    9.82      NA
ΔctDNA=POSITIVE 12     10   10.8    4.44      NA

event_summary <- circ_data %>%
  group_by(ΔctDNA) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ΔctDNA, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA Kinetics Post/Under treatment", ylab= "Overall Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("Decreased ΔctDNA", "Increased ΔctDNA"), legend.title="")

summary(KM_curve, times= c(0, 12, 24, 36))

Call: survfit(formula = surv_object ~ ΔctDNA, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ΔctDNA=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     13       0    1.000   0.000        1.000        1.000
   12      7       5    0.615   0.135        0.308        0.818
   24      3       3    0.352   0.139        0.112        0.607
   36      2       0    0.352   0.139        0.112        0.607

                ΔctDNA=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     12       0    1.000   0.000      1.00000        1.000
   12      6       6    0.500   0.144      0.20848        0.736
   24      2       3    0.222   0.128      0.04111        0.492
   36      1       1    0.111   0.101      0.00701        0.378

circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ΔctDNA, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ΔctDNA, data = circ_data)

  n= 25, number of events= 18 

                 coef exp(coef) se(coef)     z Pr(>|z|)
ΔctDNAPOSITIVE 0.6947    2.0031   0.4806 1.445    0.148

               exp(coef) exp(-coef) lower .95 upper .95
ΔctDNAPOSITIVE     2.003     0.4992    0.7809     5.138

Concordance= 0.598  (se = 0.061 )
Likelihood ratio test= 2.1  on 1 df,   p=0.1
Wald test            = 2.09  on 1 df,   p=0.1
Score (logrank) test = 2.17  on 1 df,   p=0.1

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 2 (0.78-5.14); p = 0.148"

circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels = c("NEGATIVE", "POSITIVE"), labels = c("Decreased", "Increased"))
circ_data$OS.Event <- factor(circ_data$OS.Event, levels = c("FALSE", "TRUE"), labels = c("Alive", "Deceased"))
contingency_table <- table(circ_data$ΔctDNA, circ_data$OS.Event)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 0.58792, df = 1, p-value = 0.4432

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.3783
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  0.363132 39.390335
sample estimates:
odds ratio 
  2.985191 

print(contingency_table)

           
            Alive Deceased
  Decreased     5        8
  Increased     2       10

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA kinetics post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "Vital Status",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Alive" = "blue", "Deceased" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#Association of ctDNA kinetics post/during-ICI with BOR

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ΔctDNA!="",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels = c("NEGATIVE", "POSITIVE"), labels = c("Decreased", "Increased"))
circ_data$RESIST <- factor(circ_data$RESIST, levels = c("CR", "PR", "SD", "PD"))
contingency_table <- table(circ_data$ΔctDNA, circ_data$RESIST)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test

data:  contingency_table
X-squared = 15.385, df = 3, p-value = 0.001516

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.0003154
alternative hypothesis: two.sided

print(contingency_table)

           
            CR PR SD PD
  Decreased  4  5  3  1
  Increased  0  0  3  9

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA kinetics post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "BOR",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("CR" = "blue", "PR" = "lightblue", "SD" = "orange", "PD" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#PFS by ctDNA clearance post/during-ICI

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))

survfit(Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

                  n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1  6      2     NA   19.68      NA
ctDNA.Dynamics=2 19     18   3.68    2.37    7.95

event_summary <- circ_data %>%
  group_by(ctDNA.Dynamics) %>%
  summarise(
    Total = n(),
    Events = sum(PFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="PFS - ctDNA clearance Post/Under treatment", ylab= "Progression-Free Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("Clearance", "No Clearance"), legend.title="")

summary(KM_curve, times= c(0, 12, 24, 36))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0      6       0    1.000   0.000        1.000        1.000
   12      4       1    0.833   0.152        0.273        0.975
   24      2       1    0.625   0.213        0.142        0.893
   36      1       0    0.625   0.213        0.142        0.893

                ctDNA.Dynamics=2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     19       0   1.0000  0.0000      1.00000        1.000
   12      2      17   0.1053  0.0704      0.01777        0.284
   24      1       1   0.0526  0.0512      0.00359        0.214
   36      1       0   0.0526  0.0512      0.00359        0.214

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 25, number of events= 20 

                            coef exp(coef) se(coef)     z Pr(>|z|)   
ctDNA.DynamicsNo Clearance 2.085     8.048    0.767 2.719  0.00655 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance     8.048     0.1243      1.79     36.19

Concordance= 0.672  (se = 0.049 )
Likelihood ratio test= 11.61  on 1 df,   p=7e-04
Wald test            = 7.39  on 1 df,   p=0.007
Score (logrank) test = 9.78  on 1 df,   p=0.002

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 8.05 (1.79-36.19); p = 0.007"

circ_data$PFS.Event <- factor(circ_data$PFS.Event, levels = c("FALSE", "TRUE"), labels = c("No Progression", "Progression"))
contingency_table <- table(circ_data$ctDNA.Dynamics, circ_data$PFS.Event)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 7.2505, df = 1, p-value = 0.007088

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.005477
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
    1.796359 1833.444857
sample estimates:
odds ratio 
  27.59073 

print(contingency_table)

              
               No Progression Progression
  Clearance                 4           2
  No Clearance              1          18

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA clearance post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "Progression",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("No Progression" = "blue", "Progression" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#OS by ctDNA clearance post/during-ICI

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))

survfit(Surv(time = circ_data$OS.months, event = circ_data$OS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.months, event = circ_data$OS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

                  n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1  6      2     NA   20.50      NA
ctDNA.Dynamics=2 19     16   9.82    7.75    26.7

event_summary <- circ_data %>%
  group_by(ctDNA.Dynamics) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA clearance Post/Under treatment", ylab= "Overall Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("Clearance", "No Clearance"), legend.title="")

summary(KM_curve, times= c(0, 12, 24, 36))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0      6       0      1.0   0.000        1.000        1.000
   12      5       0      1.0   0.000           NA           NA
   24      2       2      0.6   0.219        0.126        0.882
   36      1       0      0.6   0.219        0.126        0.882

                ctDNA.Dynamics=2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     19       0    1.000  0.0000       1.0000        1.000
   12      8      11    0.421  0.1133       0.2037        0.625
   24      3       4    0.189  0.0942       0.0503        0.396
   36      2       1    0.126  0.0813       0.0222        0.325

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 25, number of events= 18 

                             coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.DynamicsNo Clearance 1.6133    5.0196   0.7578 2.129   0.0333 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance      5.02     0.1992     1.137     22.17

Concordance= 0.65  (se = 0.047 )
Likelihood ratio test= 6.57  on 1 df,   p=0.01
Wald test            = 4.53  on 1 df,   p=0.03
Score (logrank) test = 5.52  on 1 df,   p=0.02

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 5.02 (1.14-22.17); p = 0.033"

circ_data$OS.Event <- factor(circ_data$OS.Event, levels = c("FALSE", "TRUE"), labels = c("Alive", "Deceased"))
contingency_table <- table(circ_data$ctDNA.Dynamics, circ_data$OS.Event)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 3.6032, df = 1, p-value = 0.05767

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.03241
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   0.8988842 150.9697807
sample estimates:
odds ratio 
   9.34674 

print(contingency_table)

              
               Alive Deceased
  Clearance        4        2
  No Clearance     3       16

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA clearance post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "Vital Status",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Alive" = "blue", "Deceased" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#Association of ctDNA clearance post/during-ICI with BOR

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$RESIST <- factor(circ_data$RESIST, levels = c("CR", "PR", "SD", "PD"))
contingency_table <- table(circ_data$ctDNA.Dynamics, circ_data$RESIST)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test

data:  contingency_table
X-squared = 14.309, df = 3, p-value = 0.002513

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.001129
alternative hypothesis: two.sided

print(contingency_table)

              
               CR PR SD PD
  Clearance     3  3  0  0
  No Clearance  1  2  6 10

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA clearance post/under-treatment", 
       x = "ctDNA", 
       y = "Patients (%)", 
       fill = "BOR",
       caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("CR" = "blue", "PR" = "lightblue", "SD" = "orange", "PD" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Progression label size

#Multivariate cox regression for PFS

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$cStage <- factor(circ_data$cStage, levels = c("0-II", "III-IV"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + cStage + CPS.Scorev2, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for PFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression for PFS v2

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional", "Metastatic"))
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative", "Positive"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + Inclusion.status + p16.status + CPS.Scorev2, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for PFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression for PFS v3

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional", "Metastatic"))
circ_data$Location <- factor(circ_data$Location, levels = c("Oropharynx/Oral", "Larynx", "Unknown"))
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative", "Positive"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + Inclusion.status + Location + p16.status + CPS.Scorev2, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for PFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Univariate PFS cox regression for variables included in MVA

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]

circ_data$ctDNA.Dynamics <- NA
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance")) #univariate for ctDNA clearance post-treatment
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data)
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 25, number of events= 20 

                            coef exp(coef) se(coef)     z Pr(>|z|)   
ctDNA.DynamicsNo Clearance 2.085     8.048    0.767 2.719  0.00655 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance     8.048     0.1243      1.79     36.19

Concordance= 0.672  (se = 0.049 )
Likelihood ratio test= 11.61  on 1 df,   p=7e-04
Wald test            = 7.39  on 1 df,   p=0.007
Score (logrank) test = 9.78  on 1 df,   p=0.002

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 8.05 (1.79-36.19); p = 0.007"

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
cox_fit <- coxph(surv_object ~ Age, data=circ_data) #univariate for age
summary(cox_fit)

Call:
coxph(formula = surv_object ~ Age, data = circ_data)

  n= 25, number of events= 20 

        coef exp(coef) se(coef)      z Pr(>|z|)
Age -0.02566   0.97467  0.02392 -1.073    0.283

    exp(coef) exp(-coef) lower .95 upper .95
Age    0.9747      1.026      0.93     1.021

Concordance= 0.512  (se = 0.079 )
Likelihood ratio test= 1.12  on 1 df,   p=0.3
Wald test            = 1.15  on 1 df,   p=0.3
Score (logrank) test = 1.16  on 1 df,   p=0.3

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 0.97 (0.93-1.02); p = 0.283"

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$cStage <- factor(circ_data$cStage, levels = c("0-II", "III-IV")) #univariate for Stage
cox_fit <- coxph(surv_object ~ cStage, data=circ_data)
summary(cox_fit)

Call:
coxph(formula = surv_object ~ cStage, data = circ_data)

  n= 25, number of events= 20 

               coef exp(coef) se(coef)    z Pr(>|z|)
cStageIII-IV 0.9241    2.5196   0.6372 1.45    0.147

             exp(coef) exp(-coef) lower .95 upper .95
cStageIII-IV      2.52     0.3969    0.7226     8.785

Concordance= 0.585  (se = 0.047 )
Likelihood ratio test= 2.55  on 1 df,   p=0.1
Wald test            = 2.1  on 1 df,   p=0.1
Score (logrank) test = 2.24  on 1 df,   p=0.1

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 2.52 (0.72-8.78); p = 0.147"

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional", "Metastatic")) #univariate for Stage/disease status
cox_fit <- coxph(surv_object ~ Inclusion.status, data=circ_data)
summary(cox_fit)

Call:
coxph(formula = surv_object ~ Inclusion.status, data = circ_data)

  n= 25, number of events= 20 

                              coef exp(coef) se(coef)      z Pr(>|z|)
Inclusion.statusMetastatic -0.1110    0.8949   0.4582 -0.242    0.809

                           exp(coef) exp(-coef) lower .95 upper .95
Inclusion.statusMetastatic    0.8949      1.117    0.3646     2.197

Concordance= 0.556  (se = 0.056 )
Likelihood ratio test= 0.06  on 1 df,   p=0.8
Wald test            = 0.06  on 1 df,   p=0.8
Score (logrank) test = 0.06  on 1 df,   p=0.8

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 0.89 (0.36-2.2); p = 0.809"

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative", "Positive")) #univariate for p16 status
cox_fit <- coxph(surv_object ~ p16.status, data=circ_data)
summary(cox_fit)

Call:
coxph(formula = surv_object ~ p16.status, data = circ_data)

  n= 24, number of events= 20 
   (1 observation deleted due to missingness)

                     coef exp(coef) se(coef)    z Pr(>|z|)
p16.statusPositive 0.2431    1.2752   0.6402 0.38    0.704

                   exp(coef) exp(-coef) lower .95 upper .95
p16.statusPositive     1.275     0.7842    0.3636     4.472

Concordance= 0.491  (se = 0.027 )
Likelihood ratio test= 0.14  on 1 df,   p=0.7
Wald test            = 0.14  on 1 df,   p=0.7
Score (logrank) test = 0.14  on 1 df,   p=0.7

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.28 (0.36-4.47); p = 0.704"

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$Location <- factor(circ_data$Location, levels = c("Oropharynx/Oral", "Larynx", "Unknown")) #univariate for disease location
cox_fit <- coxph(surv_object ~ p16.status, data=circ_data)

Warning: Loglik converged before variable  2 ; coefficient may be infinite. 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ p16.status, data = circ_data)

  n= 25, number of events= 20 

                         coef  exp(coef)   se(coef)      z Pr(>|z|)
p16.statusPositive  2.431e-01  1.275e+00  6.402e-01  0.380    0.704
p16.statusUnknown  -1.815e+01  1.309e-08  6.641e+03 -0.003    0.998

                   exp(coef) exp(-coef) lower .95 upper .95
p16.statusPositive 1.275e+00  7.842e-01    0.3636     4.472
p16.statusUnknown  1.309e-08  7.641e+07    0.0000       Inf

Concordance= 0.526  (se = 0.041 )
Likelihood ratio test= 3.49  on 2 df,   p=0.2
Wald test            = 0.14  on 2 df,   p=0.9
Score (logrank) test = 1.92  on 2 df,   p=0.4

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = -18.15 (0.78-76405588.46); p = 0.64"

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1")) #univariate for CPS score
cox_fit <- coxph(surv_object ~ CPS.Scorev2, data=circ_data)
summary(cox_fit)

Call:
coxph(formula = surv_object ~ CPS.Scorev2, data = circ_data)

  n= 25, number of events= 20 

                coef exp(coef) se(coef)     z Pr(>|z|)
CPS.Scorev2<1 0.3866    1.4720   0.4905 0.788    0.431

              exp(coef) exp(-coef) lower .95 upper .95
CPS.Scorev2<1     1.472     0.6793    0.5628      3.85

Concordance= 0.564  (se = 0.062 )
Likelihood ratio test= 0.59  on 1 df,   p=0.4
Wald test            = 0.62  on 1 df,   p=0.4
Score (logrank) test = 0.63  on 1 df,   p=0.4

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.47 (0.56-3.85); p = 0.431"

#Multivariate cox regression for OS

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$cStage <- factor(circ_data$cStage, levels = c("0-II", "III-IV"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$OS.months, event = circ_data$OS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + cStage + CPS.Scorev2, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for OS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression for OS v2

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
    ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
  ))

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional", "Metastatic"))
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative", "Positive"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + Inclusion.status + p16.status + CPS.Scorev2, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for OS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)